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Abstract 



^ Dynamic loadings produce high stress waves leading to the spallation of ductile materials 

1^3 ■ such as aluminum, copper, magnesium or tantalum. The main mechanism used herein to 

g . explain the change of the number of cavities with the stress rate is nucleation inhibition, as 

q | induced by the growth of already nucleated cavities. The dependence of the spall strength 

1 ■ and critical time with the loading rate is investigated in the framework of a probabilistic 

t-h . model. The present approach, which explains previous experimental findings on the strain- 

^ | rate dependence of the spall strength, is applied to analyze experimental data on tantalum. 

o : 
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1 Introduction 
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The impact of a projectile on a target generates two shock waves propagating in opposite directions. 
Meeting free surfaces, these shock waves reflect back as two release waves, which generally meet together 
at a definite location, the spall plane. Their superposition produces a triaxial tensile ramp loading that often 
results, prior to fracture, in the nucleation, growth, and coalescence of microvoids in most metals. This 
phenomenon is known as "ductile spalling." 
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Fig. 1. Example of a 5 mm thick tantalum sample damaged by a symmetric impact at 270 m/s by a CuC2 flyer plate 
(the shock wave traveled from top to bottom); zone I: uniaxially loaded, zone II: biaxially loaded, zone III: rapidly 
unloaded. Only the left half of the target is shown. The right edge of the image is close to the symmetry axis. 

Although discovered long ago and studied by many authors (see Meyers and Aimone, 1983; Curran et 
al., 1987; Grady, 1988 for reviews), its modeling still raises open questions. Since the pioneering works of 
Carroll and Holt (1972) and of Glennie (1972), void growth has by far been the main concern. This led many 
authors to derive elastic-viscoplastic damage models using the overall porosity as damage variable (see, 
e.g., Eftis and Nemes, 1991; Cortes 1992), often comparable to the quasi-static class of Gurson-like models 
(Gurson, 1977; Tvergaard, 1999). In these models, nucleation and coalescence are generally dealt with in 
an empirical fashion. In the recent years, however, renewed attention has been paid to these processes. The 
present paper aims at addressing the question of nucleation, coalescence being put aside for future work 
(the interested reader may refer to some recent works on this topic by Thomason, 1999; Tonks et al., 2001; 
Bontaz-Carion and Pellegrini, 2006). 

Some recent interrogations in relation to the definition of a dynamic representative volume element (Roy, 
2003; Dragon and Trumel, 2003) seem to indicate that the overall porosity is not a sufficient parameter, and 
that the entire void size distribution should be accounted for. The question of micro-inertia, neglected for 
a long time, is the subject of a continued effort (Ortiz and Molinari, 1992; Tong and Ravichandran, 1995; 
Wang and Jiang, 1997; Molinari and Mercier, 2001 ; Wu et al., 2003; Roy, 2003). Not only does it slow down 
the growth of individual voids, but also does it confine each void within an evolving neighborhood bounded 
by an elastic relaxation wave. Hence, dynamic void interactions are strongly linked to intervoid spacing, 
itself driven by the nucleation process. The latter thus appears as a crucial mechanism. This is all the more 
the case that Roy (2003), studying pure tantalum over a large range of shock levels and durations, showed 
extreme size distributions to be present in recovered samples, indicating that nucleation is a continuous 
process taking place up to coalescence. Fig. Q] shows a tantalum sample recovered after an impact at 270 
m/s by a copper flyer plate (Roy, 2003), and containing isolated voids up to about 100 fim in diameter (even 
larger voids can be observed at lower impact velocities). A detailed account of the nucleation process is 
clearly beyond the present state of knowledge, although much progress is being made using atomistic tools 
(see in particular Rudd and Belak, 2002). However, the probabilistic approach is an interesting alternative, 
as shown by Grady and Kipp (1979, 1980) and Denoual and Hild (2000) for dynamic fragmentation of 
brittle materials, and more recently by Molinari and Wright (2005) and Czarnota et al. (2006, 2008) for 
ductile spalling. In both cases, the purely deterministic description of void growth is combined with a 
stress-dependent probability of void nucleation, in the form of a Weibull-like model. Czarnota et al. (2006, 
2008) defined a probability of nucleating new voids; in addition, Denoual and Hild (2000) used a spatial 
distribution of crack nuclei among which new cracks are activated. Void interactions are also treated in a 
different fashion. Czarnota et al. (2006, 2008) used the overall porosity to describe the weakening effect of 
the already present voids, whereas Denoual and Hild (2000) considered microcrack growth as a spatially 
bounded relaxation process that inhibits nucleation in relaxed zones. In this respect, the degree of coupling 



2 



a b 




Space Space 



Fig. 2. Simplifying assumptions of the nucleation model (nucleation conditions met at arrows), (a) Physical situation: 
applied tensile stress a m and nucleation level o" nuc as continuous random fields, (b) Simplification: uniform applied 
stress, and field of nucleation thresholds split up into discrete sites of random locations and threshold values. 

is stronger in the last approach. 

It is intended here to assess the relevance of inhibition concepts for fragmentation (Mott, 1947; Grady and 
Kipp, 1980; Denoual and Hild, 2000) to analyze ductile spalling processes. Rather than precisely describing 
joint nucleation, growth processes and their couplings, this paper aims at setting the fundamentals of the 
theory to demonstrate its potentialities within the simplest possible theoretical framework. Section [2] shows 
how the deterministic and probabilistic parts of the model are interlinked, and puts the emphasis on inertial 
growth, which drives the inhibition process. Section l2.3.2l presents an application to ramp loading, generally 
agreed to be representative of the real loading in the spall plane in the lack of any phase transition process, 
and ends up with a closed-form solution of the whole problem. A very simple overall damage model is 
proposed in Section I3TT1 and yields an analytical expression for the spall stress, i.e., the maximum tensile 
stress the material can sustain during the whole spalling process. Through a thorough examination of the 
experimental data of Nicollet et al. (2001), Roy (2003), and Bontaz-Carion and Pellegrini (2006) on pure 
tantalum, the model is identified and discussed in Section HI and applied tentatively to other materials in 
Section [51 



2 Nucleation and growth in ductile materials 

2. 1 Model outline 

As introduced above, the physical process of nucleation and growth during early stages of ductile spal- 
lation is complex. Wave propagation induces transient macroscopic stress fields. At a finer (mesoscopic) 
spatial scale, local fields experience fluctuations due to the poly crystalline nature of the materials consid- 
ered (Fig. [2k). When the local hydrostatic stress cr m (x, t) exceeds some local nucleation threshold cr nuc (x), 
cavities are nucleated and start to grow. 

As shown by Roy and Villechaise (Roy, 2003), in pure tantalum nucleation sites are primarily located 
at grain boundaries, especially triple points 0] Growing cavities in turn induce relaxation zones in which 
local stresses decrease, thus decreasing the probability of nucleating voids in these zones, and out of which 

2 This picture is valid when the material is pure. When second-phase particles or precipitates are present, the so-called 
heterogeneous nucleation processes take place. The present paper focuses on the first mechanism. 
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local stresses remain unaltered. Hence, any volume element in which macroscopic stresses are uniform 
prior to nucleation evolves into a volume containing growing perturbed zones in an otherwise unperturbed 
uniformly loaded matrix. 

According to Roy (2003), isolated voids remain spherical from very small to very large sizes, implying 
that local fluctuations of material properties do not seem to influence void growth. Hence, a first simplifi- 
cation will consist in neglecting the effects of the polycrystalline nature of the material of the matrix, and 
therefore on macroscopic stresses. We thus assume uniform loading, in a pristine matrix material that con- 
tains a random spatial distribution of void nuclei at which the elastic-plastic properties of matrix material 
strongly fluctuate around their bulk value (Fig. Wp)- Furthermore, the matrix is assumed perfectly plastic 
hereafter. Neglecting temperature, viscosity and strain hardening is performed for the sake of simplicity, 
and can be relaxed in more detailed (future) analyses. During further evolutions, activated voids are the 
only local heterogeneities that will affect macroscopic stresses. In this context, a voided volume is viewed 
as a matrix loaded by a uniform hydrostatic tensile stress cr m , containing several (possibly overlapping) 
perturbed zones. 

Second, the joint effects of local stress fluctuations and of local weaknesses are accounted for through a 
stress-dependent nucleation probability. It will be further assumed that the inhibition phenomenon is total 
in strongly relaxed zones. Hence, matrix stresses will be considered as the only driving force for nucleation 
and growth. 

Third, given the high level of triaxiality, as well as the spherical shape of the voids observed by Roy 
(2003), macroscopic shear stresses will be neglected, such that = a m 5ij. From now on, a m will simply 
be referred to as "the stress." 

Growth drives the extension of relaxation zones, and thus the inhibition process. The growth model must 
therefore be carefully chosen. On the one hand, as stressed by Ortiz and Molinari (1992), Wang and Jiang 
(1997), Roy (2003), Dragon and Trumel (2003), Molinari and Wright (2005), and Czarnota et al. (2006), 
inertial effects are overwhelmingly important. On the other hand, elasticity should not be neglected, since it 
has a strong effect on early growth (Denoual and Diani, 2002; Roy, 2003). Advantage will be taken here of 
a simplified approach proposed by Roy (2003), from the work of Forrestal and Luk (1998). This approach 
shows that growth cannot take place if the macroscopic stress is less than a cavitation threshold, as shown 
by many authors in the quasi-static case (Mandel, 1966; Hou and Abeyaratne, 1992; Denoual and Diani, 
2002). 

We now proceed to assemble the above-listed ingredients. In a pristine examination volume V subjected 
to uniform stress a m (t), we assume the number N of active nucleation sites, of associated random nucleation 
stress cr nuc (x) where x is the site location, to follow a point-Poisson distribution of intensity n tot (the average 
volume density of active sites). The probability of finding iV active sites in V is 

P(N,V)= { ^^-exp(-n tot V). (1) 

In the above definition, a nucleation site at location x is said active at t (i.e., can potentially nucleate a 
void) if a m (r) > a nuc (x) for any past time < r < t. It will effectively give birth to a void only if not 
inhibited (effects of inhibition are dealt with in Section [2.3. II ). Introduce then cr max (t) = max < r < i cr m (r), 
the maximum hydrostatic stress reached up to time t. According to experimental findings (Roy, 2003), the 
density of nucleated cavities is stress-dependent. This prompts us to further write P(N, V) in the form of 
the so-called Weibull-Poisson law by taking (Gulino and Phoenix, 1991; Jeulin, 1991; Denoual and Hild, 
2002) 
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"totw = ™o , (2) 

where m is the Weibull modulus which characterizes the scatter in nucleation levels (weak scatter corre- 
sponds to a high m value, and conversely), cr is a scale parameter relative to a reference density n , and 
(*) are Macauley brackets that denote the positive part of *. In Eq. (OQ), the product n tot V thus represents 
the average number of sites in V where cr m has overcome the nucleation threshold. Eq. © indicates that the 
higher (T max (t), the more nucleation sites are active. It should be noted that a classical Weibull expression is 
retrieved within the weakest link framework, see Appendix lAl 

Since Eqs. (OQ) and © describe the probability of activating N sites in a pristine uniformly loaded volume 
V, they also hold (with V replaced by V) in the uniformly loaded part V of a larger voided volume, by 
definition of a m , and by the above assumption of total inhibition. The volume V is found by subtracting 
from V the volume of inhibited zones, thus accounting for possible overlaps between individual inhibition 
zones that grow out of each activated site. Since inhibition is related to stress relaxation, V depends on the 
growth model, which is addressed now. 

2.2 A simplified growth model 



a b 
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Fig. 3. Equivalent hollow sphere model, (a) Real elasto-plastic hollow sphere, (b) Simplified representation of equiv- 
alent elastic energy, (c) Schematic representation of overlapping relaxed zones. See Section l3~Tl for a discussion of (b) 
and (c). 

Cavity nucleation can be understood as a bifurcation process in the sense of Hou and Abeyaratne (1992). 
Once nucleated, any new cavity starts to grow. As shown, for example, by Hopkins (1960), Hunter and 
Crozier (1968), Glennie (1972), or Roy (2003), an isolated growing cavity of radius a(t) can be seen as 
an expanding volume bounded by an elastic relaxation wave at radius b(t). This volume consists in an 
outer elastic zone, and an inner elastic-plastic region, separated by an evolving boundary of ("plastic") 
radius c(t) (Fig. [3k). Both regions are referred to as "the matrix" hereafter. Denoual and Diani (2002) and 
Tonks et al. (2001) showed that the early growth can be decomposed into three distinct phases. The first 
one is essentially elastic, until the hydrostatic stress reaches a "cavitation threshold" (see below). There, 
bulk elastic energy release induces a violent elastic-plastic expansion of the cavity, until the third phase of 
stationary expansion is established. 

Strong relaxation occurs inside the elastic-plastic zone. This is illustrated numerically, by submitting 
a hollow tantalum sphere of initial outer and inner radii of 12.5 and 0.58 fim, respectively (i.e., an initial 
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Fig. 4. Space-time diagram for a hollow sphere (initial outer radius: b = 12.5 /jm, initial inner radius a = 0.58 pm) 
submitted to a stress rate of 2 GPa/^s on the outer radius. The material parameters are those of tantalum (TableQ]). The 
plastic zone coincides with the region in which the local pressure decreases, in spite of an overall pressure increase, 
(a) Plastic zone evolution and (b) stress rate indicator. 

porosity of 10~ 4 ) to a hydrostatic stress ramp of 2 GPa ps^ 1 applied on the outer boundary. Fig. |4] shows 
the space-time domain where the matrix is yielding (Fig. 4a), and that of varying a m (Fig. 4b). It is seen 
that a m decreases inside the plastic zone although the applied stress keeps increasing. Plastic zones can thus 
be seen as (and identified to) inhibition zones for further void nucleation, and this is exploited in the next 
section. This shows that unlike previous works (Wu et al., 2003; Molinari and Wright, 2005; Czarnota et al., 
2006), it is not sufficient to establish a link between the macroscopic stress and the cavity radius a, namely 
the link between these quantities and c must be known as well. 

Roy et al. (2002) and Roy (2003) checked numerically that the cavitation stress is independent of the 
macroscopic strain rate, and that the transient regime is brief. Accordingly, and since this allows for closed- 
form solutions, a purely stationary model is used here, with a(t) = at, c(t) = ct, where a and c are constant 
growth velocities. The approach used (detailed in Appendix |B]) is adapted from the work of Forrestal and 
Luk (1988), itself derived from earlier works dealing with isolated cavity growth under internal hydrostatic 
stress (Hopkins, 1960; Hunter and Crozier, 1968). 

Thus, for an isolated cavity in an infinite medium subjected to a remote tensile stress cr m (t), an implicit 
equation links c and a m (t) to a (Eqs. (|B.14I) and (IB.15I) ). A first-order expansion valid in the low stress rate 
regime (assuming c <C c P and a <ti cp, where cp - 



' Kj p is the so-called plastic velocity (Zel'dovich 
and Raizer, 2002), K is the bulk modulus, and p is the reference density), then provides the additional 
proportionality relationship 



-1/3, 



(3) 



where (3 is defined by Eq. (IB. 16b in terms of K, p the shear modulus, and Y the yield stress. For most 
materials [3 1. In turn, a similar first-order expansion provides relationship (IB. 1 81) . namely, 



Oo(o"m/O"cav ~ 1) 1/2 



(4) 



between the void growth velocity and the applied tensile stress, where a cav is the cavitation threshold, and 
where d is a characteristic void growth velocity in the material. Both quantities depend only on K, p and Y, 
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with d depending on p as well (see Appendix |B] for explicit expressions of these quantities). Combining 
Eqs. © and © yields 



c = kc P (cr m /cr cav - 1) 1/2 , (5) 

where k = /5~ 1//3 d /cp is a numerical coefficient (k c P being a characteristic growth velocity of the plastic 
region). For Al, Cu and Ta, k varies between 0.3 and 0.5. 

Eq. ([5]) constitutes a particular instance of the more general class of threshold-like expressions 

c = kc P (a m /a cav - l) a , (6) 

where a > is a stress- sensitivity exponent, and where the nucleation stress a nuc of Fig. [2] is identified 
to the cavitation threshold o" cav . No significant growth of the micro void population should occur unless 
cavitation conditions are met. This general expression covers the present case, as well as the "quasi-brittle" 
case for which a = (Denoual et al., 1997). In the case of monotonically increasing loading cr m (i), upon 
integrating © over time we obtain c(t) in the form 

c(t)=C(t-t mc ), (t>W), (7) 
where C is some function and t nuc the nucleation time obtained as a solution to 

C"m(^nuc) "'cav (8) 

2.3 Elementary cell assembly 

2.3.1 Dynamic inhibition model 

So far, we described the behavior of isolated cavities only, in a deterministic way. The collective behavior 
of the population of voids is now considered. Henceforth, overlined quantities are used for macroscopic 
variables that represent statistical (or more phenomenological) averages of their microscopic counterparts. 

The intrinsic probabilistic nature of the nucleation and growth process should be embodied in some 
random variability of the local elastic and plastic properties of the material Y, //, K and p . Eq. © shows 
that under some prescribed time-dependent loading, a random set {cr cav } of cavitation or generic nucleation 
thresholds (see Fig. |2]) can be mapped to a random set {t nuc } of nucleation times. Randomness in the 
process is thus introduced through the following crucial bold assumption that emphasizes the part played 
by nucleation times, namely, in Eq. © the nucleation time t nuc , which physically depends on the above 
material parameters and on the local loading, will be considered as a random variable, whereas material 
parameters, and parameters that define the field loading function, will be considered as "averaged" ones 
whenever they enter the definition of the function C itself. 

Section [2T21 substantiates the identification between plastic regions and zones of total nucleation inhibi- 
tion. Accordingly, the inhibition volume V inh associated to an isolated cavity is taken hereafter proportional 
to the plastic radius c to the third power 

V inh = V inh (t - t imc ) =Sc 3 , (9) 

where S is a shape parameter, and the functional time dependence of V inh stems from Eq. ©. 
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Fig. 5. Inhibition and horizon concepts, (a) Inhibition phenomena. In grey are relaxed zones where void nucleation is 
inhibited, (b) Horizon of a site P. Any active site in the grey zone inhibits further cavity nucleation at P. 

New voids will nucleate from active nucleation sites (in the sense of Sec. 12. II) only if they do not belong 
to any relaxed zone produced by previously nucleated growing voids, as depicted in Fig. 5a. Thus n nuc , the 
volume density of the centers of nucleated voids, is related to n tot defined in Eq. © by 



dt 



Pi 



dn t , 



ot 



inhJ 



dt 



(10) 



with n nuc (0) = n tot (0) = 0. This equation, which implements inhibition effects in the model, involves the 
inhibition probability (identified with an overall volume fraction of inhibited regions) 



-Pinh(^) 



exp 



-V inh (t)n tot {cr m (t)} 



where V^nh, the mean volume of the inhibition zone, is defined by 



V- mh {t)n tot {o m (t)} 



V inh (t - t)^^ {o- m (r)} dr. 
dr 



(ID 



(12) 



Eqs. CCD) and (|T2|) (Denoual et al., 1997), which originate from the Poisson hypothesis, Eq. (OQ), are derived 
in AppendixO which makes clear that Eq. CCD) accounts for the overlaps of inhibition zones (the derivation 
uses the horizon concept described in Fig. 5b, which constitutes another way to look at the inhibition 
process). From the point of view of mathematical morphology, this model constitutes an instance of a 
Boolean islands model (Jeulin and Jeulin, 1981; Serra, 1982). Also, in the context of isothermal diffusive 
phase transformations, the three latter equations are central to the Kolmogorov-Johnson-Mehl-Avrami 
(KJMA) kinetic theory of nucleation and growth (Kolmogorov, 1937; Johnson and Mehl, 1939; Avrami, 
1941). Eqs. CCD and (fT2l are valid for any density n tot and any shape of interaction zones of volume 
Vi n h. The present framework is thus adaptive to incorporate many different inhibition phenomenologies. In 
particular, the same approach can be used to analyze dynamic fragmentation of brittle materials (Grady and 
Kipp, 1979, 1980; Denoual and Hild, 2000, 2002). In that case, inhibition is induced by stress relaxation 
around propagating cracks, as was also the case for the shell fragmentation problem studied by Mott (1947). 

Bearing in mind the particular time dependence of V[ n h in Eq. ©, it is observed that the time-integration 
in Eq. (fT2l) is over the nucleation time. According to our above hypothesis of considering the nucleation 
time as a random variable, Eq. (fT2l) indicates that its associated probability density at time t imposed by the 
Weibull-Poisson process © is (with r > 0) 
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V(t nuc = r;t) = - tttt-^ Wm(r)} , (13) 

That Wm{t)} dr 

where 9 is the Heaviside step function. 

Finally, an expression for the average void volume fraction / in the examination volume is obtained as 
follows. Eq. (0) implies, via a = /5 1//3 c, the following proportionality relationship between the individual 
cavity volume V C!W oc a 3 and the corresponding inhibition volume V? n h oc c 3 

V cw = PV iBh . (14) 

Using Eq. (fl2]>. the average cavity volume K cav follows as 

F cav (t)=/5F inh (t). (15) 

Since the individual voids and inhibition zones are of same centers, they obey the same statistics. The 
porosity / is thus 

fit) = 1 - exp [-V cav (t)n tot {a m (t)}] (16) 
and simply relates to the inhibition probability by 

/ = 1 - (1 - Pinhf- (17) 

This relationship is illustrated by Fig. 3c, interpreting in the present context white zones as voids of overall 
volume fraction /, and dotted zones as inhibited zones of overall volume fraction P inh . 

2.3.2 Application to ramp loading 

In general the number of nucleated cavities must be computed numerically. The nucleation Eq. (fTOl) 
involves the matrix stress in the non-inhibited zones, a m . The link with the overall stress is given in Sec- 
tion 13.11 The computation is particularly simple for the particular case of ramp-stress loading a m = at 
with constant stress-rate a that yields a closed-form solution of practical interest for experimental analyses. 
Upon integrating Eq. © over time, and introducing the nucleation time i nuc = o c&N ja according to the first 
paragraph of Section l2.3.1l the individual inhibition volume © reads 



Mnh — S 



kcp ( a 



a + 1 \a. 



(t - t 



nuc ) 



(18) 



for t > t nuc , and zero otherwise. The corresponding cavity volume V CSM follows from Eq. (fT4l . 

At this stage, it proves useful to introduce a dimensionless flaw density h = n/n c , time t = t/t c , volume 
V = V/V c and stress cr m = a m /a c . Two ways of defining those dimensionless quantities are relevant here. 
Both are based on the condition 

n c V c = 1 (t = t c ) (19) 

that expresses the fact that some characteristic volume V c contains on average one site at time t c . 

Computing P ilin requires identifying V c with the inhibition volume, whereby the above condition reads 

n C iVci = 1, n ci = n tot [a m (t ci )}, V ci = V inh (t ci ), (20) 
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where the subscript ci denotes characteristic quantities associated to inhibition. A characteristic stress is 
defined by a ci = at ci . From Eqs. and (1201) . the characteristic parameters follow as 



( n _l_ 1 ^3_m_.3a 1 l/[m+3(a+l)] 
I"" 1 " 1J O °cav 



n (kc P ) 3 Sa m+3a 
(a + l)nr i)/m ^av^ 



3m/[m+3(a+l)] 
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(« + 1) V "V C 3 « a 3 



l/[m+3(a+l)] 

n (fccp) 3 5" 

Upon carrying out the integration in Eq. (fT2)) . Eq. (fTT|) reads 



(21) 



rm+3(a+l) 



-Pinh = 1 — exp — B(m, 3(a + l))t 



where B is a modified Euler function of the first kind 

Jo T(p + q + l) 

and the closed-form solution of Eq. (flOl) yields 

-m/[m+3(o+l)] 



(22) 



(23) 
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m B(m, 3(a + 1) 



m + 3(a + 1) 
where 7 is the incomplete gamma function 
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(24) 



7(p,x) 



(25) 



so that 7(p, a; — > +00) = T(p). Eq. ([24]) is the eracf solution to Mott's problem (1947) extended to three- 
dimensional cases with an initial flaw density modeled by a power law function. Fig. [6] shows the change of 
the dimensionless density 77 nuc with the dimensionless time t. At early times t < 1, virtually no inhibition 
is observed, i.e., P in h « and 77 nuc « 77 tot . Conversely, at late times 1, P in h « 1 and saturation occurs. 
The higher the Weibull modulus m, the higher the density at saturation (Fig. [7]). 

Computing / instead requires identifying V c with the void volume, whereby the characteristic parameters 
obey 



^■cc^cc 1) ^cc ^tot [^m(^cc)] > ^cc ^cav^cc)) 



(26) 



where the subscript "cc" denotes characteristic quantities associated to cavities. Similarly, the characteristic 
stress is defined by a cc = crt cc . Then, Eqs. (fT4l) and (|26|) provide 



, o-l/[m+3(a+l)] T/ = ta. W[m+3(a+l)] 

CI 1^ *) CC Cl/ 1 ^ * 



-l/[m+3(a+l)] 



(27) 



and the overall porosity /, Eq. (fT6l) . takes on the form 
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Fig. 6. Dimensionless nucleated density n nuc , Eq. (|24l) . vs. dimensionless time t for three different Weibull moduli m 
when a = 1/2 (ra c denotes either n c ; or n cc ). 
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Fig. 7. Dimensionless nucleated density at saturation, n nuc (oo), vs. modulus m when a = 1/2 (n c denotes either n c 
or n cc ). 
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B(m,3(a + l))i 
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(28) 



Remark that Eq. dTTl) follows from replacing k by /c/3 1//3 in Eq. (f2TT) . The only quantitative difference be- 
tween Eqs. (T28l) and (T22l) resides in the definition of the characteristic parameters (i.e., t = t/t ci for inhibition 
and t = t/t cc for cavities). These results are exploited below in the framework of a simplified constitutive 
model. 
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3 A simplified constitutive model 



3.1 Homogenization approach for dynamic loadings 

In usual homogenization approaches to the computation of the overall constitutive law of disordered 
porous media, some void spatial distribution is prescribed in advance, all voids being by hypothesis in 
mutual long-range elastic interaction, and the homogenization problem amounts to finding suitable ap- 
proximation schemes for these interactions. Such approaches quite generally provide estimates of stress 
fluctuations in the matrix (due to pore elastic interactions), which can be considered as evenly spread in 
the latter. In stark contrast with this situation, the dynamical impact conditions considered here consist in 
loading a pristine matrix with a uniform stress state <j m in the first place, this initial state being perturbed 
afterwards by relaxation waves originating from nucleated growing voids. As a consequence, stress fluctu- 
ations in the matrix are more localized (at least until significant overall relaxation is achieved through some 
"percolation" of the relaxed zones), and it should be clear that standard homogenization techniques ought 
not be straightforwardly transposed to this case. 

The following alternative two-step approach is adopted instead, motivated by the elastic decoupling of 
the voids in the first stages of the spall process. In a first step, the elementary voided elastic-plastic cell of 
radius b, with void radius a, in which the stress is heterogeneous but equal to a m on its boundary (Fig. 3a), 
is replaced by an equivalent cell of radius b containing a fictitious void of radius a cq (region of null stress), 
outside which the stress is uniform and equal to cr m (Fig. 3b). The volume fraction of fictitious void in the 
equivalent cell being written S(c/b) 3 , where 5 is an unknown proportionality constant, it is proposed here 
to compute 5 by requiring the elastic energy densities in the real and fictitious systems to be equal. The 
equation for 5 thus reads 



1 



2 



la 



2 



-(a : IT 1 : <7) cell = 1 - 5{c/bY (29) 



2 K 



where (-) ce ii denotes a volume average over the elementary cell and where C is the usual (isotropic) tensor 
of elastic moduli built on K and fi. The l.h.s. of Eq. (1291 ), which involves microscopic hydrostatic and shear 
stress components, can be computed using the stress of the exact solution for the real elastic-plastic cell, 
derived in AppendixEJ Since this solution also provides a and b in terms of a m , the outcome is an expression 
of 5 as a function of a m . The associated fictitious void volume V eq = Sa 3 q = 5V in h is then introduced (it is 
recalled that V inh = Sc 3 ). 

Since fictitious voids obey the same point-Poissonian statistics as real ones, the same token that was used 
to relate / to Pi n h given the relationship between V cav and Vi n h in Sec. 12.3. 1 L can be re-used here to relate 
P in h to an overall fictitious porosity / eq given the above relationship between V| n h and V C!W . This second 
step provides the macroscopic relationship analogous to Eq. (fTTT) 

1 - f eq = (1 - P inh ) 5 . (30) 

In the macroscopic equivalent system, the stress outside the fictitious voids is now homogeneous every- 
where, equal to cr m (Fig. |3fc). Hence the expression of the macroscopic stress a m in terms of a m reads 

Om = (1 - feq)(T m = (1 ~ Pinh)V m , (31) 

where P in h and 5 depend on a m . In this relation / eq plays the part of an overall damage variable D in the 
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standard relation a m = (1 — D)a m of damage theory (Lemaitre and Chaboche, 1990). To emphasize this 
connection the notation D = f eq is used from now on. 

The computation of 5 from the solution of Appendix |B] is quite involved. Besides, Eq. (|29l ) is not free 
from arbitrariness since other energetic equivalences could be proposed that explicitly involve an additional 
kinetic energy term as proposed by Wang and Jiang (1994) and Molinari and Mercier (2001). Bearing in 
mind the present exploratory purpose, a pragmatic and simplified approach is preferred that consists in 
considering only the limiting cases 5=1, whereby D = P in h, and 5 = (5 whereby D = f. These limits, 
respectively, provide upper and lower "pseudo-bounds" (if not rigorous ones) to D. The former assumes 
that relaxation is total in the plastic zone, and neglects elastic relaxation, such that the equivalent volume 
is the plastic zone volume. The latter neglects any relaxation, such that the equivalent volume is the void 
volume. The relevance of these "bounds" is established below in Sec. 14.31 by comparison to experimental 
results. 

For ramp loading, the constants required to write down in dimensionless form the equations of type 
(1311) that stem from each "bound" have been worked out in Sec. 12.3.21 An example of the dimensionless 
macroscopic stress a m / a c as a function of the dimensionless time t, which reduces to the same master curve 
for both "bounds", is displayed in Fig.[8j 




Fig. 8. Dimensionless macroscopic stress provided by Eq. (T3TT) . vs. dimensionless time when m = 8 and a = 1/2 {a c 
denotes either a c \ or u cc ). 

Completing the above approach in order to arrive at a full constitutive relationship between a m and the 
macroscopic strain e m requires making additional assumptions, and is not needed here. 

3.2 Spall criterion and spall strength 

The spall strength is the quantity of primary interest in dynamic ductile damage experiments. It is defined 
as the maximum macroscopic stress 7x s sustained by the material during the damage process. Given the 
relationship a m = a m (a m ) between the macroscopic stress <r m and the microscopic stress a m in non- 
perturbed, uniformly loaded, regions of the matrix (see Section I2T21) . the macroscopic spall strength a s can 
be obtained as a s = a m (a s ), where a s is the microscopic spall stress solution of 
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p±(a m = a B )=0. (32) 
da 



in 



The spall strength a s corresponds to the maximum stress in the plot of Fig. [8l For the ramp load solution, 
Eqs. (1221) and (1311 ), the derivative in (1321 is carried out using t = cr m /a and 5 = 1 or [3 in the solution. It 
vanishes for a dimensionless critical time for spall 

t s = {[m + 3(a + l)]B{m,3(a + l))y 1/lm+3(a+1)] (33) 
and a corresponding macroscopic spall strength 

a s = a c {[m + 3(a + l)]B(m, 3(q + l))e}~ 1/Im+3(oH " 1)] , (34) 

where e = exp(l), £ s = t/t cc , a c = a cc for the upper "bound", and t s = t/t ci , a c = a ci for the lower 
"bound". At the spall point, the damage parameter is equal for both "bounds" to 

D s = l-exp{-l/[m + 3(a + l)]}. (35) 



4 Analyses of experiments on tantalum 

4. 1 The material 

Tantalum is a transition metal of great interest for studying dynamic ductile damage mainly because 
of its high mass density (16,660 kg/m 3 ), good dynamic strength and very high ductility in wide strain 
rate and temperature ranges. The samples used herein are machined from 5 mm thick cross-rolled and 
fully recrystallized plates. Advanced elaboration process and heat treatment resulted in a very high purity 
material (99.98 wt%). The main (embrittling) impurities are 15 wt ppm O, 15 wt ppm C and less than 10 
ppm N, with a homogeneous micro structure characterized by equiaxed grains of typical size 90/im, and a 
weak residual texture. Either optical microscopy, SEM or SIMS examinations did not reveal any localized 
heterogeneity down to a ~ 5 /xm scale, namely no second-phase hard particle nor impurity gradient at 
grain boundaries. The lack of preferable nucleation sites has been revealed by dynamic tensile tests on 
smooth and notched axisymmetric samples, where failure does occur in any case by ultimate thinning of 
the elongated ligament rather than through inclusion-induced damage, for stress triaxialities ranging from 
0.3 to 1 (Roy, 2003). This material is consequently an almost ideal polycrystal for studying homogeneous 
ductile nucleation (Roy, 2003). 

Mechanical properties of tantalum have been carefully determined from ultrasonic measurements, quasi- 
static and dynamic uniaxial testing on both as received and shocked material (Roy, 2003). During the release 
stage following the initial shock compression, tantalum behaves roughly as an isotropic elastic perfectly 
plastic medium (Juanicotena, 1998; Roy, 2003). This holds both at the macroscopic scale during release 
wave interaction when no damage occurs and at the mesoscopic scale around growing voids, where high 
strain-rate gradients are roughly balanced by thermal softening at large strain. The relevant properties of 
tantalum in the range of stress and strain states of interest are summarized in TableQ] 

Twenty-two plate impact experiments (Nicollet et al., 2001; Roy, 2003; Llorca and Roy, 2003; Bontaz- 
Carion and Pellegrini, 2006) were performed and/or analyzed for the present paper. Impact velocity, flyer 
plate material and flyer plate thickness were selected as relevant parameters for varying both shock pressure 
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Parameter 


Value 


Compressibility modulus K 


191 CPa 


Shear modulus ft 


69GPa 


Mass density p B 


16,660kg/m 3 


Yield stress V 


700 MPa 



Table 1 

Tantalum material parameters. 



and pulse duration, and are summarized in Table 2. This essentially induces variations in the position of 
the plane of maximum tensile stress (the spall plane) and in the mean and maximum achievable tensile 
stress state along this plane. The diagnostics used to study the condition for damage and spall to occur are 
Doppler laser interferometry to record the velocity of the target free surface (overall structural response of 
the sample plate) and qualitative and quantitative metallurgical analyses of the soft-recovered samples. The 
most significant results derived from this microstructural examination have been reported elsewhere (Roy, 
2003; Llorca and Roy, 2003; Nicollet et al., 2001; Bontaz-Carion and Pellegrini, 2006). 



300 




0.5 1 1.5 2 2.5 3 3.5 4 
Times ([is) 



Fig. 9. Example of a free-surface velocity record exhibiting a pull-back signal (Roy, 2003). 
4.2 Data extraction 

Time-resolved in situ measurements in the spall plane are not yet possible, and an inverse methodology 
must be adopted. As a result, data extraction is performed from free surface velocity records. In the spall 
plane, progressive damage induces local relaxation waves whose macroscopic consequence is a so-called 
pullback signal (see Fig. [9]). For first-order estimations of relevant data (spall plane location, spall strength, 
critical time to fracture), a simple analytical elastic method is often used (Romanchenko and Stepanov, 
1980). This method proves successful at low shock pressure (lower than the material dynamic yield strength) 
or at high pressure (when elastic behavior can be neglected regarding plastic hydrodynamic component) 
(Meyers, 1994). This is definitely not the case for tantalum, whose dynamic yield strength (or Hugoniot 
Elastic Limit) is known to be less than an order of magnitude lower than its spall strength in the range of 
loading paths of interest. 
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Accurate data extraction requires an analysis of the complex wave pattern induced by the plate impact. 
One- and two-dimensional numerical simulations have consequently been performed using the Lagrangian 
explicit hydrocode Hesione (a proprietary code of the Commissariat a l'Energie Atomique). In order to 
extract the thermomechanical fields in the region of interest (the spall plane) as accurately as possible, a 
tabulated equation of state and a Preston et al. (2003) viscoplastic constitutive law are used (Juanicotena, 
1998). 

Within the relatively low shock pressure range investigated here (low temperature increase and weak 
plastic strain during shock and release at the macroscopic scale), these relationships do predict an essen- 
tially elastic perfectly plastic behavior during unloading, consistent with the analytical parameters summa- 
rized in Table \T\ These relationships fitted from dedicated experimental databases on shock and uniaxial 
compression behavior of this tantalum grade (Roy, 2003), yield very good correlation with experimental 
results used in this study regarding shock and release behavior. A fracture criterion is added, leading to 
instantaneous mesh opening at a given tensile stress threshold (spall strength), and fitted numerically for 
each simulated experiment. 

This numerical procedure is sufficient to extract the following data from free-surface velocity records: 
the stress rate in the matrix, the critical time and the spall strength. In order for the extraction procedure to be 
as accurate as possible, two features are particularly sought in matching numerical results and free-surface 
velocity records, namely the minimum velocity preceding pullback signajjO and the subsequent ringing 
velocity frequency, which suggests efficient prediction of both spall plane position (which was compared 
with the experimental value for some experiments), effective maximum tensile stress and associated critical 
time, in a far more accurate way than using the simplified analytical method presented by Roy (2003). In 
this fashion, critical time (spall criterion activation) and mean tensile stress rate are derived from numerical 
stress history prediction at the spall plane before fracture. The corresponding values are summarized in Ta- 
ble |2l Associated error estimations are derived from numerical investigation of impact velocity uncertainty, 
mesh size, and artificial viscosity sensitivity of the calculated spall strength for some typical experiments. 
For some of these experiments, quantitative relevant damage activation measures were also derived using 
metallurgical observation of sample slices coupled with optical profilometry and image analysis (Roy, 2003) 
for an estimation of the three-dimensional damage state. In particular, the volume density of nucleated voids 
could be measured in the vicinity of the spall plane. These values are also given in Table |2] 

4.3 Identification and validation 

Fig. [10] shows the change of volume density of nucleated voids n tot , Eq. ©, as a function of the shock 
pressure. This plot is restricted to data obtained from shots A1-A5, A7, A9, A12 and B8 only. These shots 
involve only moderate pressures so that void coalescence presumably remains limited. Moreover, in the 
impact configurations considered, the shock pressure is equal to the negative of a m , the maximum stress 
in the matrix, which takes place in non-inhibited regions that exist whenever coalescence is marginal. The 
assumption of a constant stress rate pulse (i.e., ramp load) is applied to tantalum to determine the Weibull 
parameters of Eq. ©. The best power-law fit displayed in Fig. [Jj)]provides an exponent m = 8, a moderate 



3 We emphasize that a completely fractured plane at the macroscopic scale is not a necessary condition for pullback- 
type free surface velocity evolution, as highlighted by Llorca and Roy (2003) and Roy (2003). Primary internal 
energy release leading to pullback velocity (early re-acceleration) has been experimentally shown to be initiated in 
the vicinity of the spall plane at a given low incipient damage level. This is consistent with the basic hypothesis of the 
spall criterion developed in part 3. 
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Shot Experimental 



Impact simulations 





Target 


Imp. nature 


Imp. thickn. 


Impact vel. 


Void density 


Shock press. 


Spall strength 


Stress rate 


Critical time 




thickn. 




(mm) 


(ms- 1 ) 


(number/ mm 3 } 


(CPa) 


(GPa) 


(GPa/ps) 


(MS] 




(mm) 




















±0.01 




±0.01 


±2.0% 


±30% 


±5% 


±5% 


±12% 


±17% 


Al 


4.95 


CuC2 


3.00 


252 


256 


5.5 


N/A 


N/A 


N/A 


A2 


4.95 


CuC2 


1.00 


271 


60 


5.8 


N/A 


N/A 


N/A 


A3 


4.95 


CuC2 


2.00 


270 


N/A 


5.87 


5.32 


17 


0.36 


A4 


4.95 


CuC2 


2.00 


270 


75 


5.87 


N/A 


N/A 


N/A 


A5 


4.95 


CuC2 


4.00 


268 


301 


5.88 


5.37 


11 


0.48 


A6 


4.95 


CuC2 


3.00 


269 


N/A 


5.9 


5.32 


15 


0.34 


A7 


4.95 


Ta 


3.00 


207 


340 


5.93 


5.1 


15 


0.32 


AS 


4.95 


Ta 


3.00 


303 


N/A 


8.9 


6.6 


42.8 


0.15 


A3 


4.95 


Ta 


3.00 


306 


3660 


9.01 


6.6 


49.7 


0.12 


A10 


4.95 


Ta 


3,00 


307 


N/A 


9.01 


6.6 


53.6 


0.12 


All 


4.95 


Ta 


3.00 


415 


N/A 


12.5 


7.5 


81.1 


0.08 


A12 


495 


Ta 


3.00 


424 


150,000 


12.7 


N/A 


N/A 


N/A 


Bl 


3.96 


AU4G (Al 2017) 


3.97 


412 


N/A 


4.92 


N/A 


9.1 


N/A 


B2 


3.92 


AU4G (Al 2017) 


3.49 


481 


N/A 


5.86 


5.35 


17.3 


0.38 


B3 


3.96 


AU4G (Al 2017) 


1.99 


570 


N/A 


6.5 


5.9 


16.2 


0.34 


B4 


3.94 


AU4G (Al 2017) 


3.01 


551 


N/A 


6.89 


6.15 


26,9 


0.24 


B5 


3.93 


AU4G (Al 2017) 


3.50 


654 


N/A 


8.3 


7.05 


44.7 


0.16 


B6 


3.96 


AU4G (Al 2017) 


4.01 


528 


N/A 


8.55 


7.1 


47 


0.15 


B7 


3.99 


AU4G (Al 2017) 


3.99 


671 


N/A 


9.06 


7.2 


106.4 


0.Q7 


BS 


3.95 


AU4G (Al 2017) 


Q.49 


1039 


2250 


9.06 


7.2 


106.4 


0.07 


B9 


3.90 


AU4G (Al 2017) 


3.49 


900 


N/A 


11.95 


7.9 


99.9 


0.07 


BID 


3.96 


AU4G (Al 2017) 


3.4B 


1110 


N/A 


15.25 


7.05 


80 


0.04 



Table 2 

Parameters of shock experiments. Shots A1-A12 (resp. B1-B10) are those of Roy (2003) and Llorca and Roy (2003) 
(resp. Nicollet et al., 2001; Bontaz-Carion and Pellegrini, 2006). Third row: standard uncertainties. Sixth column: 
voids densities measured by image analysis on recovered samples. N/A indicates unavailable data. 

value indicative of weak scatter in nucleation levels. A value of a = 700 MPa for the scaling stress equal 
to elastic limit is used, whence the density n = 7.9 x 10~ 6 mm -3 is obtained. Upper and lower theoretical 




1 10 100 



Shock pressure (GPa) 

Fig. 10. Volume density of pores ntot vs. shock pressure for tantalum. The solid symbols are experimental points and 
the line is the best fit of Eq. (f2]). 

"bounds" for the critical time vs. stress rate obtained from ([33]) with a — 1/2 and m = 8 are displayed in 
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u 

1 




0.01 



Stress rate (GPa/iis) 



Fig. 11. Critical time vs. stress rate. The solid symbols are experimental points and the dashed lines are the "bounds" 
built from Eq. (T33T ) with a = 1/2, and m = 8 determined from Fig.flOl 





Stress rate (GPa/Lis) 



Fig. 12. Spall strength vs. stress rate for tantalum. The solid symbols are experimental points and the dashed lines are 
the bounds built from Eqs. (I33T ) and (l34l ) with a = 1/2, and m = 8 determined from Fig.flOl 



Fig. QT] Almost all experimental points are seen to lie within these "bounds". Besides the overall trend is 
consistent with the slope of the latter. Likewise, upper and lower theoretical bounds derived from (|34l ) with 
a = 1/2 and m = 8 are compared to experimental data in Fig. [121 in log-log scale, which illustrates the 
power-law increase of the spall strength with the stress rate. Though the experimental points are linearly 
correlated with a slope lower than that of the bounds, they lie between the latter in the considered range of 
loadings, which is quite satisfactory. Thus, the rate sensitivity of the spall strength can be described by the 
present model with no need to incorporate a time-dependent constitutive equation of the matrix. 
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5 Analyses of data on aluminum and magnesium 

Kanel et al. (1996) performed experiments on aluminum and magnesium. In both cases, the spall strength 
was shown to be approximated by a power-law of the strain rate. In the present analysis, as in all the devel- 
opments derived herein, the effect of the temperature is ignored. Consequently, only experiments performed 
at ambient are considered. By using a s = a t s = t s a t c oc a t c , any of Eq. (f2TT) or (|27l) for t c vs. &, and the 
proportionality a oc e (of elastic origin, and legitimate in non-relaxed regions of uniform <r m ), the following 
strain-rate dependence is obtained for the microscopic spall strength 

a s oc eP (36) 

with 

3 

V = — r^T — \ 77; (3V) 

where e denotes the average strain rate in the experiments. In this expression, the only unknown is the mod- 
ulus m, provided a value of a = 1/2 is chosen as in the previous experiments on tantalum. For aluminum, 
a value 77 = 0.059 is found, which would lead to a value of m = 46 and for magnesium, 77 = 0.072 so that 
m = 37. These two (high) values of m are an indication of a small scatter in terms of nucleation level when 
compared to tantalum (Table [3]) for which a gradual and scattered nucleation was observed. 



Parameter 


Tantalum 


Aluminum 


Magnesium 


a 


0.5 


0.5 


0.5 


Weibull modulus m 


8 


46 


37 



Table 3 

Nucleation parameters for tantalum, aluminum and magnesium. 



6 Conclusion and perspectives 

We proposed a probabilistic model for nucleation and growth in ductile fracture, using Poisson-Weibull 
statistical concepts, which are usually applied to brittle materials. We showed, through analyses of several 
sets of experimental data in spall experiments, that these concepts are well suited to describing ductile 
fracture as well. In particular, we arrived at a simple explanation for the power-law dependence of the 
spall strength vs. strain rate observed by Kanel et al. (1996). The proposed model makes use of a velocity- 
dependent extension of the concept of cavitation stress in metals. Though it has been presented, for sim- 
plicity, in the framework of ideal plasticity, expressions of cavitation thresholds that account for hardening 
are available (Mandel, 1966; Bishop and Hill, 1945), and could be easily appealed to. Investigations of the 
influence of hardening on the present nucleation theory, as well as that of viscoplastic behavior, are left to 
future work. Also, next steps should consist in implementing the full model in a finite-element hydrocode, 
and in extending its range of validity to the coalescence regime. 
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Appendix 
A Weibull distribution 

The probability of finding at least one nucleation site (i.e., the "weakest link") in a uniformly loaded 
domain ft is 

P{N > 1, ft) = 1 - P{N = 0, ft) = 1 - e - Vn °(W m )/a r^ (A.l) 
When the domain is not uniformly loaded, we instead have 

P(N > l,ft) = 1 - e - v ^M(°M)/vor ; ( A .2) 
where V e g denotes the effective volume (Davies, 1973) 

with <jm = max<j m (x). (A. 3) 

Eqs. (IA.1I) and (IA.2I) are the Weibull model (Weibull, 1951) written in the context of ductile damage (see 
also Czarnota et al., 2006). 

B Derivation of Eqs. © and ® 

B. 1 Preliminaries 

We detail here the steps leading to Eqs. © and ©, in the stationary growth regime studied by Forrestal 
and Luk (1988). In this one-dimensional spherical approach, a cavity of radius a(t) grows at constant 
velocity d in an infinite elasto-plastic medium submitted to an initially uniform hydrostatic stress state 
cr m (t). This growth perturbs the stress field within a partially relaxed volume of radius r = b(t) = c^t, 

where cl = \J (K + 4/i/3)/p is the velocity of longitudinal elastic waves. The front b(t) separates the outer 
medium at rest in a state of uniform stress, from the inner perturbed region expanding with the growing void. 
The inner region is divided into an external elastic shell c(t) < r < b(t), and a shell at yield that surrounds 
the cavity, a(t) < r < c(t). 

In the steady-state growth regime where c(t) — ct, a self-similar solution for the radial displacement u 
is sought for in the form u(r, t) = c(t) There, £(r, t) = r/c(t) is the scaled radial coordinate, and w(£) 
is the scaled displacement. Moderate stress is assumed so as to neglect (i) density variations in the elastic 



shell, (ii) non-linear elasticity, and (iii) convection terms (Forrestal and Luk, 1988). Using £ = — £ c/c, the 
velocity and acceleration read 

ii = (u- £u')c, (B.l) 
u = i 2 u"c 2 jc. (B.2) 



eff 



/ d 3 a; 
Jq 



0"M 
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Eq. (IB. 1 1) provides the scaled velocity = u{r, t)/c. 

For further use, we introduce the scaled yield stress, shear modulus, applied hydrostatic stress, and 
density, respectively, as y = K/if, g = 2/i/K, a(t) = a m (t)/K, p(£) = p(r, t)/p , where p is the 
reference material density. In usual metals, 

y « 5 < i. (B.3) 

5.2 Elastic shell 

Combining linear elasticity relationships and the momentum equation 

d r a r + (2/r)(ov - a e ) = pii, 

where oy and erg are, respectively, the radial and hoop stresses, and introducing ^ L = c/cl, yields the 
differential equation 

i - ii e) u " + (2/e) cee'-s) =o. (b.4) 



Its solution is of the form «(£) = C'- 1 ^ + C^l — 37 2 £ 2 )/£ 2 , where the integration constants are 
found from boundary conditions. The first one is w(£ = — &m/ {3Kjl)> an d stems from the applied 

external boundary traction. The second one is ~~ ^'(Ok=i = ^/ (2/J-), which expresses the yield 

condition ag — a r = Y (tensile case) at the elastic-plastic boundary. The solution for £ G (1, l/7z) is then 



a F (1- 7 lQ 2 (1 + 2 7 l 
3 6 + 6a* (1-71K 2 



Denoting by v the Poisson ratio, the corresponding radial stress reads 

~e s 2y (1 - 7l{) [(1 - 2v)(\ + 7L fl + (1 + tQTggj 

= ff " T (i-2„)(i- 7 ye ' <B ' 6) 

B.3 Plastic shell 

Mass conservation, namely, d t p + [d r + (2/r)](up) = 0, provides 

v' + (2/Z)v=(t-v)p'/p. (B.7) 



Introduce now the plastic velocity cp = y K/p , and (after j L ) another scaling of c as 7 P = c/cp. The 
yield condition o~g — a r = Y, combined with linear elasticity in the form Tr a = a r + 2ag = 2>K (p / P — 1), 
gives 

d r a r = -K (p /p 2 ) d rP . (B.8) 
Using (IB.8I) in the momentum equation then provides 

ip'=(liev'p-2y)f. (B.9) 
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Eliminating p'(£) from Eqs. (IB.7I) and (IB .91) then yields 

1 " e 7p P 2 (f - «0l u ' + 2 (1 - y p) (v/O = -2yp. 



(B.10) 



Eqs. (lB.9l ) and (IB. 101) constitute a system for p and t>, should variations in p be accounted for. Upon neglect- 
ing their higher order influence in v at moderate stress, and assuming y <C 1, see (|B.3I) . Eq. (|B.10I) reduces 
to 



(l-T* P t 3 )v' + 2{v/t) 



■2y. 



(B.ll) 



Note that this equation also assumes that v(£) <C £, which is satisfied if the material velocity is much 
lower than the velocity of the void boundary. However, finite-element calculations of void expansion (Roy, 
2003) indicate that this assumption is expected to hold everywhere except near the void boundary where the 
velocity gradient is highest. The difference induced by neglecting this term on the overall behavior is small 
anyway (Forrestal and Luk, 1988; Roy, 2003), see Fig. lB.ll below. 

Continuity of the material velocity at the elastic-plastic interface provides the boundary condition u(l) = 
y/g. Then, the solution of Eq. (IB .111) in the interval £ 6 (a/c, 1) is 



v p (0 



y 



lie 



2 £2 



1-7?. 



;i+7?>/<?)-£ 



+ 



y 



;i-7i>r)iog 



;i + 7 pQ(i-7p) 
;i- 7 p0(i + 7p)' 



(B.12) 



With a P = d/cp the scaled void growth velocity, the radial stress in the same interval reads, upon integrat- 
ing (IB .81) and using (IB.7I) under the above approximations 
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5.4 Complete and approximate solutions 



The equations for the void growth velocity then consist in the relations ap = 7pf p (£ = a/c) and 
5^(1) = cr^(l). The first one reads 



gfc/y + 1 7p/g + l , 1, 1 + ap 1-7P 

2~ = 7P-; 2" + o lo § 71 77777 7' 

1-ap I-Tp 2 (l-a P )(l + 7 p) 

whereas, setting k = c^/cp = (1 + 2g/3) 1 / 2 , the second one yields 



(B.14) 
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n + ± log ii±Ml^ +log ^-<* 
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ap(l 



7p) 



. (B.15) 



Seeking low-order expansions of Eqs. (|B.14I) and (IB. 151) . a is computed as a function of a m by first 
looking for a solution of Eq. (IB. 141) in the perturbative form 7p = J7fc>i A k a P , where the unknowns A k are 
determined order-by-order. To leading order in a P , the solution is 
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Fig. B.l. Dimensionless velocity 7p = c/cp of the plastic zone vs. dimensionless growth velocity up = a/cp of a 
cavity in a compressible and incompressible medium (the value of cp at finite compressibility is used for all curves). 
Material parameters of tantalum (Table [D, except for the bulk modulus in the incompressible case. 

7 ^ r vs ap , ^rtLtmy + n (B , 6) 

9{y + 3/2) 2/i 3K 

where the approximated value of (3 stems from (|B.3I) . Next, inserting the expansion into Eq. (IB. 151) . assum- 
ing a relationship a(ap) = a c + Y,k>i Bk<yp where a c and the Bk are unknowns, and again simplifying the 
coefficients with (IB.3I) . yields 

or = ^(1 - log/3) + [2 - O ((^ 2 ) 1/3 )] o? P + O (4) , (B.17) 

where the orders of the neglected terms are indicated. Hence B 2 = 2 in the incompressible limit. We do not 
reproduce its full expression, quite involved but easily retrieved with a symbolic calculator. The first term 
in the r.h.s. is the scaled cavitation stress, a cav = a cav /K, first computed by Bishop et al. (1945)0 and later 
on by Mandel (1966) for finite compressibility under the form a cav = (2Y/3) {1 + log £/ [3(1 - v)Y}}, E 
being Young's modulus. 

Growth occurs only if a m > cr cav . Hence from (IB. 1 VI ). for a m > a cav , the pore growth velocity behaves 

as 



a ~ 



<ioK/<T cav - 1) 1/2 , (B.18) 



where do = [er C av/(-E?2 po)] 1 ^ 2 is a characteristic pore growth velocity of the material with B 2 — 2. Using 
the full expressions of (3 and B 2 , we obtain for Al, Cu and Ta: d — 289, 224 and 145 m/s respectively, 
and cr cav ~ 0.11, 0.89, and 2.75 GPa respectively. For comparison purposes, we note that cp ~ 5092, 3589, 
3386 m/s for these materials, respectively, so that d is lower than c P by more than one order of magnitude. 
Neglecting compressibility provides, with B 2 = 2, values of d lower than the above ones by a relative error 
of about 5 x 10~ 3 . Though it is strictly valid for a constant applied stress (since a = const, by hypothesis), 



However, their expression, written in terms of \x and Y, is that of the incompressible case. 
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Eq. (IB. 1 81) nonetheless provides the leading-order behavior for a time-varying stress, which is the type of 
loading considered in Sec. 12.21 where use is made of this equation. In practice, transient corrections mainly 
consist in damped oscillations around this leading behavior, as was checked by finite-element calculations, 
and are neglected in this work. 

In the incompressible limit (where cl = cp = oo) Eq. (IB. 14b reduces to a = [Y/ (2/i)] 1 / 3 c (this relation 
is encapsulated in the equations of Carroll and Holt (1972) in the limit of vanishing initial porosity). Hence 
in general, it is expected that a <C c < cp. Combined with the above low-velocity solution, this suggests 
the following approximation to (IB. 14b 

a P = /3 lp /(l- lp ), (B.19) 

which contains in particular the incompressible limit (where 7p — > 0). This approximation, which preserves 
(IB. 1 Vb up to the neglected terms, and which can be solved analytically for jp, proves useful to compute 
numerically the stress in numerical implementations of the model. 

Fig. IB . 1 1 compares Eq. (IB. 14b with either stress- or velocity-driven finite-element numerical results. 
These data points are reasonably well reproduced by the solution of Eq. (IB. 14b . in spite of the underlying 
approximations. The solution to Eq. (IB. 1 9b is indistinguishable to the eye from the latter. Also shown are 
the linear approximation (IB. 1 6b and the above incompressible (linear) solution. 

It should be noted that the incompressible limiting value B 2 = 2 markedly differs from the value 
B 2 = 3/2 which one easily deduces from Carroll and Holt's (1972) incompressible calculation in the 
limit of zero initial porosity, where convection is accounted for. Though a detailed study of the influence of 
convective terms in the compressible case lies beyond the scope of this paper, this difference indicates that 
convection may be important in accurately determining the coefficient d in (IB. 1 8b . the difference between 
the approaches concerning a numerical coefficient of order one. Taking B 2 = 3/2 instead of 2, Eq. ( IB. 17b is 
compatible with the work of Molinari and Wright (2005) in the limit of stationary growth of incompressible 
materials, and close to the result given by Tonks et al. (2001). Thus, the obtained cavitation threshold and 
the general form of this law hold in any case, which is a sufficient conclusion for the present purpose. 



C Inhibition probability 

To define the probability that a point x at a time t be relaxed, it is preferable to invert the problem by 
looking into the past of the considered site to know if a cavity is able to inhibit its nucleation (this method, 
first proposed by Cahn (1996), was found independently by two of the present authors (Denoual et al., 1997; 
Denoual, 1998)). Two zones are distinguished. First, a zone in which the nucleated cavities never inhibit 
the considered site (see dashed part of Fig. 5(b) when r < t). In the second (complementary) zone, any 
nucleated cavity will inhibit x. This zone is referred to as the horizon (Cahn, 1996; Denoual et al., 1997; 
Denoual, 1998). 

The inhibition probability Pmh(t) is written as the product of the elementary probabilities AP^(r) 
1-Pinh(t)= IlAP ? (r), (Cl) 

T=0 

where AP^(r) is the probability of finding no new sites during a time increment At in a zone VJ n h(t — r). 
It suffices to apply Eq. (OQ) with V = V[ n h(£ — r) for an intensity (dn tot /dr){o"(r)}Ar, since it still is a 
Poisson process 
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AP^(t) = exp 



dn 



'tot 



{a(T)}ArV inh (t - t) 



dr 



(C2) 



The probability Pmhit) becomes 



1 - Pinh(t) = exp 



dn 



tot 



Mr)}Ar^(i-r) 



(C3) 



L r=o 



dr 



In the continuous limit At — > 0, rewriting the sum as an integral eventually yields Eqs. CCD) and (fl~2|) . 
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